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A macroscopically uniform model of a two-dimensional electron system is proposed to study 
nonequilibrium properties of electrical conduction. By molecular dynamics simulation, the steady 
state distribution function P y of electron velocity in a direction perpendicular to an external driv- 
ing force is calculated. An explicit form of P y is determined within the accuracy of the numerical 
simulation, which fits the numerical data well even in the regime where a local equilibrium descrip- 
tion is not valid. Although the entire structure of P y is different from that of a local equilibrium 
distribution function, the asymptotic structure of the tails of P y in the limit of large absolute values 
of the velocity is identical to that of a Maxwell distribution function with a temperature which is 
different from that in the equilibrium state and the kinetic temperature in the steady state. 

PACS numbers: 05.60.-k, 71.10.-w, 05.70.Ln 
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One of the central roles of statistical mechanics is to 
provide a principle for calculating the probability dis- 
tribution p of a microscopic state without solving equa- 
tions of motion. In equilibrium states it is given as the 
canonical distribution. In nonequilibrium steady states 
(NESSs), by contrast, such a principle is not known. 

In studies of nonequilibrium systems with some con- 
crete models, it is generally easier to compute a single- 
particle distribution function, such as a velocity distri- 
bution function (VDF), than to calculate p itself. In 
equilibrium states, the VDF is, of course, the Maxwell 
distribution. Although the VDF contains less informa- 
tion than p since the VDF is a marginal probability dis- 
tribution, it yields a restriction on p and knowledge of 
the VDF is thus expected to provide a clue to finding a 
principle of nonequilibrium statistical mechanics. 

VDFs have been investigated mainly in models de- 
scribed by kinetic equations, such as the Boltzmann 
equation and the Bhatnagar-Gross-Krook equation 
Such kinetic models are said to be valid for dilute gas 
systems. Furthermore, VDFs in most studies on those 
models are given as a power series of degrees of nonequi- 
librium (such as a temperature gradient). From these 
facts, it may be interesting to investigate a model other 
than kinetic models using a method without power series 
expansions. 

One of the approaches in such studies is molecular dy- 
namics (MD) simulation on models which obey micro- 
scopic equations of motion. VDFs have been investigated 
with MD simulations in shear-flowing systems [2j and 
thermal-conducting lattice systems 0,111] • I* 1 those stud- 
ies, deviations from local equilibrium distribution func- 
tions (LEDFs) have mainly been considered. However, 
sufficient characterization of VDFs from the viewpoint 
of statistical mechanics has not yet been achieved. An- 
other example of MD simulations was a study of the dis- 
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tribution function of energy currents carried by a sin- 
gle particle in thermal-conducting particle systems [f|. 
In that work, the authors found that the energy cur- 
rent distribution function has tails at large absolute val- 
ues of the current, which are well described by equilib- 
rium distributions with appropriate temperatures. This 
tail structure might be a universal property of steady 
state distribution functions. That is, for a wide range 
of nonequilibrium systems, including thermal-conducting 
and other systems, tails of distributions such as VDFs 
and energy current distributions might asymptotically 
approach certain equilibrium distributions, although the 
complete structures would be different from those of the 
equilibrium distributions. 

In the present paper, to demonstrate the above idea, 
we investigate a VDF in an electrically conducting sys- 
tem. We introduce a model in which a macroscopically 
uniform NESS is realized, and perform MD simulation to 
calculate a VDF in a direction perpendicular to an exter- 
nal driving force. We find an explicit form of the VDF 
within simulation accuracy. Introduction of a velocity- 
dependent "reciprocal temperature" helps to determine 
the form. 

Model. The system of the simulation is composed of 
electrons, phonons, and impurities, which are modeled as 
classical particles. A schematic diagram of the system is 
shown in the inset of Fig. [TJ The linear dimension in the 
a-direction of the system is denoted by L a (a — x, y, z). 
Electrons (whose mass, charge, and total number are 
denoted by m e , e, and N e ) can move only in a two- 
dimensional (2D) x-y plane (shown in dark gray) which 
is located on the top (z = L z ) of the three-dimensional 
(3D) system. Phonons (whose mass and total number 
are denoted by m p and N p ) can move throughout the 
entire 3D system. Impurities (whose total number is de- 
noted by Ni) are fixed at random positions, and play 
the role of a random potential. Nf B of the impuri- 
ties are uniformly distributed throughout the 2D elec- 
tron system and A^ 3D throughout the 3D system under 
the 2D system (N? D + N? D = Ni). Periodic boundary 
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FIG. 1: Average of electron velocity in the ^-direction, plot- 
ted against E. The dashed line represents the slope of the 
linear response. Inset: schematic diagram of the simulation 
system. Electrons are confined in a 2D plane (dark gray). 
Under the 2D electron system there is a three-dimensional 
phonon system, the bottom of which is a thermal wall (light 
gray). Impurities are distributed throughout all regions of 
the system. An external electric field E is applied in the re- 
direction. 



conditions are imposed in the x- and y-directions. The 
boundary conditions in the ^-direction, applied only to 
phonons, are an elastic potential wall for the top bound- 
ary (z — L z ) and a thermal wall (shown in light gray) 
with temperature To for the bottom boundary (z = 0). 
We apply an external electric field E (acting on elec- 
trons) uniformly in the x-direction to drive the system 
to a noncquilibrium state. We assume that interactions 
among all kinds of particles are present. The interaction 
potential between the fc-th and Z-th particles is given by 
<Pki = y(vnax{0, dki}) 5 / 2 . Here, y is a constant of inter- 
action strength, and dki = Rk + Ri — I'M is the overlap 
of the potential ranges. Rk is the radius of the potential 
range (R e , R P , and Ri for an electron, phonon, and im- 
purity, respectively), and r\a is the inter-center distance 
between the particles. The energy supplied from E to 
electrons is transferred to phonons by electron-phonon 
interactions and dissipates through the thermal walls, by 
which the system retains its energy balance. In the sim- 
ulations we take e, i? e , m e , a reference energy, and the 
Boltzmann constant as the units, and fix the following 
parameters: L x — L y — 200, L z 
m p = l, T = 1, N e = 1000, N p = 4000, N'f 
Nf B = 1000, and y = 4000. 

This model is an extension of the two-dimensional 
model for the MD simulation of electrical conduction, 
which has been previously proposed in Ref. [6|. A pos- 
sible experimental situation corresponding to the model 
in Ref. occurs when the edges of a 2D electron sys- 
tem are in contact with large sample holders (which have 
large heat capacity and good heat conduction). On the 
other hand, the present model is a reflection of a typi- 
cal experimental setup of a quasi-2D electron system in 
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a semiconductor device. That is, a 2D electron system 
is realized around the top of a bulk substrate (which is 
modeled by the 3D system of phonons and impurities) 
and the substrate is mounted on a large sample holder 
(which is modeled by the thermal wall). The present 
model is macroscopically uniform in the y-direction as 
well as in the x-direction, whereas the previous one was 
uniform only in the direction parallel to E. This enables 
us to investigate bulk properties in both the x- and y- 
directions. It should also be noted that this model is 
isotropic in the y-direction and the net currents of parti- 
cles and energy in this direction are absent on average. 

In Fig.Q]we show Tv-dependence of the average electron 
velocity (v x )e in the x-direction. Here (■■■)e denotes 
the average in the steady state for E. A linear response 
is observed in the regime of small E and a nonlinear 
response in the regime of large E. 

Deviation from LEDF. Because of the macroscopic 
spatial uniformity of the model, steady states in the 
model are also uniform. That is, local quantities such 
as local number densities and local kinetic temperatures 
in steady states are almost independent of the position 
in the 2D electron plane. We may therefore use electrons 
in the entire region of the plane to calculate a VDF of 
electrons. In this paper we investigate distribution func- 
tions P y {vy) of electron velocity v y in the y-direction. 
[The corresponding results in the x-direction (parallel to 
E) will be presented elsewhere 0] ; these have more com- 
plicated structures and properties than those of P y pre- 
sented here.] In Fig. [2] we show semi-logarithmic plots 
of P y , plotted against v y 2 . The data for positive v y al- 
most completely overlap with those for negative v y in 
NESSs (E = 0.1 and 0.4), as well as in the equilibrium 
state (E = 0). This indicates that P y (v y ) is symmetric 
at v y — 0, and is consistent with (v v )e — 0, even in the 
nonlinear response regime. These results are due to the 
isotropy in the y-direction of the model. 

As an LEDF in the y-direction, we employ a naive 
distribution function defined by 



fy C ( v v) = \Jrne/2irT2 exp (-m e ^/2TJ') , (1) 

where T* = m e (v v 2 )E is a kinetic temperature in this di- 
rection (this is simply a Maxwell distribution with tem- 
perature Tj 7 ). In Fig. fj] we also plot f y oc (in this plot 
zero-mean Gaussian distributions are drawn as straight 
lines). We observe that the tails of P y and f y oc are dif- 
ferent for large values of E, whereas they are almost the 
same for E = 0. Moreover in the inset of Fig. [2J we 
plot the deviation P y (v y ) — f y ° c (v y ). From this figure we 
see that the difference around v y = is also large in the 
NESSs. The symmetric behavior of the deviation is again 
consistent with the isotropic nature of the model. 

Explicit form of P y . To explore the functional form of 
P y we next investigate a Uy-dependent "reciprocal tem- 
perature", [3y(Vy), which is defined by 



(iy{Vy) 



2\n[P v (v y )/Py(0)} 



(2) 
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FIG. 2: (Color online) A semi-logarithmic plot of the distri- 
bution P y , plotted against v y 2 . The results for E = 0, 0.1, 
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and 0.4 are shown. The data for v y > 
overlapped. The dashed lines represent the LEDFs f y 
solid lines depict fitting curves of Eq. Q for E = 
0.4. Inset: deviation P y (v y ) - f y ° c {v y ) from the LEDF. Also 



shown are the fitting curves of Eq. 
for E — 0.1 and 0.4 (solid curves). 
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If P y is a Maxwell distribution, (3 y is independent of v y 
and is equal to the reciprocal temperature of the Maxwell 
distribution. In Fig. [3] we show u^-dependence of /3 y . 
Figure [3] (a) is a result for the equilibrium state (E — 
0). In this case f3 y is almost independent of v y and is 
nearly equal to the reciprocal temperature 1/Tq (= 1) 
of the thermal wall within the precision of the numerical 
simulation. This is consistent with the fact that in an 
equilibrium state, P y is a Maxwell distribution with To. 
Although the data fluctuations are rather large for i^'s 
close to zero and in the tails, they are simply due to the 
restricted accuracy of the numerical simulation in these 
regions and they can be expected to become smaller as 
we average over greater numbers of simulation samples. 

Figures (b) and (c) show the results for NESSs 
(E = 0.1 and 0.4, respectively). As \v y \ becomes large, (3 y 
varies (increasing in (b) and decreasing in (c)) monoton- 
ically and tends to converge to a certain constant value 
which is different from both 1/T and I /TV. This 
implies that P y approaches a Maxwell distribution as 
My | —> oo. Furthermore, we find that the overall be- 
havior of f3 y is well described by a Gaussian function. 
That is, (3 y is well fitted by 



Py( v v) = -Bexp(- 



-Dv y 2 )+p™, 



(3) 



where /3£°, B, and D are fitting parameters. In the insets 
of Figs. M (b) and (c), to see this more clearly, we show 
semi- logarithmic plots of \/3 y — /3£°| versus v y 2 . We ob- 
serve exponential decay behavior of \/3 y — /3£°| as functions 
of v y 2 , which indicates that the variation of \(3 y — / 
almost equivalent to that of Gaussian functions of v y . 
Although the data for v v 's close to zero and in the tails 
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FIG. 3: (Color online) Mains: (3 y for (a) E = 0, (b) E = 0.1, 
and (c) E = 0.4, plotted against v y . The dotted and dashed 
lines represent 1/T e " and /3£°, respectively (l/TJ' ~ fJ^ ~ 
1/To = 1 when E = 0). The solid lines are fitting curves of 
Eq. @. Insets: semi- logarithmic plots of (b) /3^° — f3 y and 
(c) P y — /3^°, respectively, plotted against v y 2 . The data for 
v y > and v y < are overlapped. The solid straight lines 
represent fitting curves of Eq. © (subtracted from (3^ in (b) 
and reduced by in (c), respectively). 



fluctuate rather violently, they would converge to a single 
curve as the number of samples increases. From Eqs. © 
and ^ we thus obtain an explicit functional form of P y 
as 



p y( v v) = ^(0)exp 



-Dv 



m e v y ' 



(4) 



This is a main result of the present paper. In Fig. [5] we 
plot curves of Eq. Q for E — 0.1 and 0.4 as solid lines. 
These also reveal that Eq. (j4|) is a good fitting function 
of P y for almost all values of v y . 

E- dependence. Figure 2] shows the independences of 
parameters in Eq. ((4]) . In the top of Fig. [4] we compare 
is /3^° and 1/TjF. Although ^-dependences of these two 
quantities are similar, their detailed values are different. 
fl^ > l/TJ' at small values of E (as can be seen also 
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FIG. 4: Top: ^-dependence of /3£° (circles) and 1/T* (trian- 
gles). Bottom: _B-dependence of B. 



in Fig. [3(b)), whereas < 1/T* at large values of E 
(as is seen also in Fig. 13(c)). In the bottom of Fig. [4] 
we show the result for B. There exists a value E c of the 
electric field such that B < if < E < E c and B > 
if E > E c (E c ~ 0.15 in the parameter values of the 
present simulation). This is consistent with the result 
for which is mentioned above. That is, > 1/T» 
if < E < E c and < 1/T| if £ > E c . 

Summary and discussion. In this paper we have in- 
troduced a model of a two dimensional classical electron 
system which is macroscopically uniform. In the model 
with an external driving force (electric field E), a spa- 
tially uniform nonequilibrium steady state (NESS) is re- 
alized. By using molecular dynamics simulations we have 
investigated the distribution function P y (v y ) of an elec- 
tron velocity v y in a direction perpendicular to E in a 
steady state. When the driving force is large, the local 
equilibrium (LE) description is not valid. We have deter- 
mined an explicit form of the velocity distribution func- 



tion (VDF) as Eq. Q , which is valid up to numerical ac- 
curacy, even for NESSs where the LE description breaks 
down. The functional form has a velocity-dependent 
"reciprocal temperature" /3 v (v v ). In the limit of large 
absolute values of v y , this reciprocal temperature con- 
verges to a constant value which is different from 
both the equilibrium temperature and the kinetic tem- 
perature, and the VDF thus approaches asymptotically 
a Maxwell distribution with temperature 1 //3^° ■ We have 
also found that there exists a crossing value of the driving 
force below (above) which /3£° > /3 y (0) < A,(0)). It 
might be difficult to obtain such a nontrivial dependence 
of (3 y on E by a naive perturbation expansion in terms 
of E. 

The asymptotic property of the tails of the VDFs in 
NESSs is similar to that of nonequilibrium distribution 
functions of microscopic energy currents in thermal con- 
duction models [5]. In those models the distribution 
function of energy currents parallel to the temperature 
gradient has a tail for large negative (positive) values of 
the current which asymptotically obeys an equilibrium 
distribution with temperature T_ (T+). Here, T_ and 
T + are different from the local temperature at the place 
where the distribution is measured, but equal to the local 
kinetic temperatures at slightly forward and backward 
regions, respectively. In the present model of electrical 
conduction, the tails of the VDF asymptotically obey 
an equilibrium distribution with temperature 1 / . Al- 
though l/Py™ is equivalent in the tails for positive and 
negative values of velocity (because of the isotropy of 
the model in the direction perpendicular to the driving 
force), it differs from both the equilibrium temperature 
and the kinetic temperature. We therefore expect that it 
is a universal feature in NESSs that the distribution func- 
tion of a microscopic current approaches an equilibrium 
distribution in the limit of large currents. This might 
be related to a representation of the nonequilibrium dis- 
tribution of microscopic states, which has been derived 
recently 

Unlike the results for thermal conduction models, we 
do not yet know the physical meaning of the temperature 
Also we do not understand the meanings of other 
parameters in the functional form of the VDF and the 
crossing value of the driving force. These subjects remain 
for future research. 

The author is grateful to A. Shimizu for helpful dis- 
cussions and comments. 
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